function ray_invert_ps(varargin)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function ray_invert_ps
%
% This function is a front-end to the raytrace3d PSinvert routine to jointly
% invert for P-wave velocity and P:S velocity ratio using arrivals from a
% css3.0 database. This function requires Antelope to be installed on the
% system.
%
% This function requires the presence of two other functions:
%   ray_make_traveltime.m
%   ray_latlon2xyz.m
%
% Usage:
% ray_invert_ps(database,model_file,P_pars_file,ratio_pars_file,damping,relative_weight,[subset],[locations])
%
% Required Inputs:
%   database:        The path to a css3.0 database containing the following
%                    tables: site, affiliation, arrival, assoc, origin, and
%                    event
%
%   model_file:      The starting 3D model file used in any raytrace3d routine
%
%   P_pars_file:     The path to the input_mps.txt file for P waves required
%                    by raytrace3d PS_invert (see raytrace3d help for file
%                    description)
%
%   ratio_pars_file: The path to the input_mps.txt file for P:S ratio
%                    required by raytrace3d PS_invert (see raytrace3d help
%                    for file description)
%
%   damping:         The value of the damping parameter to be used in the
%                    inversion
%
%   relative_weight: The relative weighting between P-wave velocity changes
%                    and changes to the P:S velocity ratios.  A value of 1
%                    favors both equally; values towards 0 favors changes
%                    to P and values greater than 1 favor changes to the
%                    P:S ratio
%
% Optional Input:
%   subset:          A valid string for subsetting a Datascope database
%
%   locations:       The path to the output file containing updated
%                    earthquake hypocenters generated by the ray_locate 
%                    function
%
% Output (to filesystem):
%   inversion_model.model:      The model result from the inversion
%
%   inversion_model.model.hits: The number of ray hits at each model node
%
%
% Author:
% Matt Gardine
% May 2010
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Checks that Antelope is installed on the system.
if (exist('dbopen') ~= 3)
    error('Error: Antelope must be installed on the system to use this function')
end

% Checks for the existence of the ray_make_traveltimes function 
if (exist('ray_make_traveltime') ~= 2)
    error('Error: This function is dependent on ray_make_traveltime. Please add this function into the path')
end

% Checks for the existence of the ray_defaults file
if exist('ray_defaults','file')==2
    [ref_lat,ref_lon,ref_alt,projection]=ray_defaults();
    disp('ray_defaults file found.')
    disp(['ref_lat = ' num2str(ref_lat)])
    disp(['ref_lon = ' num2str(ref_lon)])
    disp(['ref_alt = ' num2str(ref_alt)])
    disp(['projection = ' projection])
else
    ref_lat=17.01;
    ref_lon=-105.99;
    ref_alt=0;
    projection='flat';
    disp('ray_defaults file NOT found. Values used:')
    disp(['ref_lat = ' num2str(ref_lat)])
    disp(['ref_lon = ' num2str(ref_lon)])
    disp(['ref_alt = ' num2str(ref_alt)])
    disp(['projection = ' projection])
end

% Checks for the existence of ray_latlon2xyz
if (exist('ray_latlon2xyz') ~= 2)
    error('Error: This function is dependent on ray_latlon2xyz.  Please add this function into the path')
end

switch nargin
    case 6
        database=varargin{1};
        model_file=varargin{2};
        P_pars_file=varargin{3};
        ratio_pars_file=varargin{4};
        damping=varargin{5};
        relative_weight=varargin{6};
        
        db = dbopen(database,'r');
        db_site = dblookup(db,'','site','','');
        db_affil = dblookup(db,'','affiliation','','');
        db_arr = dblookup(db,'','arrival','','');
        db_assoc = dblookup(db,'','assoc','','');
        db_orig = dblookup(db,'','origin','','');
        db_event = dblookup(db,'','event','','');

        db = dbjoin(db_orig,db_event);
        
        db = dbsubset(db,'orid==prefor');
        
        db = dbjoin(db,db_assoc);
        db = dbjoin(db,db_arr);
        db = dbjoin(db,db_affil);
        db = dbjoin(db,db_site);

        db_P = dbsubset(db,'phase=~/P/');
        db_S = dbsubset(db,'phase=~/S/');
        
        ray_make_traveltime(db_P);
        system('mv traveltimes.tsv P_traveltimes.tsv');
        
        ray_make_traveltime(db_S);
        system('mv traveltimes.tsv S_traveltimes.tsv');
        
    case 7
        database=varargin{1};
        model_file=varargin{2};
        P_pars_file=varargin{3};
        ratio_pars_file=varargin{4};
        damping=varargin{5};
        relative_weight=varargin{6};
        
        db = dbopen(database,'r');
        db_site = dblookup(db,'','site','','');
        db_affil = dblookup(db,'','affiliation','','');
        db_arr = dblookup(db,'','arrival','','');
        db_assoc = dblookup(db,'','assoc','','');
        db_orig = dblookup(db,'','origin','','');
        db_event = dblookup(db,'','event','','');

        db = dbjoin(db_orig,db_event);
        
        db = dbsubset(db,'orid==prefor');
        
        db = dbjoin(db,db_assoc);
        db = dbjoin(db,db_arr);
        db = dbjoin(db,db_affil);
        db = dbjoin(db,db_site);

        db_P = dbsubset(db,'phase=~/P/');
        db_S = dbsubset(db,'phase=~/S/');
        
        if exist(varargin{7},'file')~=2
            subset=varargin{7};
            db_P = dbsubset(db_P,subset);
            db_S = dbsubset(db_S,subset);
            
            ray_make_traveltime(db_P);
            system('mv traveltimes.tsv P_traveltimes.tsv');
            
            ray_make_traveltime(db_S);
            system('mv traveltimes.tsv S_traveltimes.tsv');
        else
            locations=varargin{7};
            
            ray_make_traveltime(db_P,locations);
            system('mv traveltimes.tsv P_traveltimes.tsv');
            
            ray_make_traveltime(db_S);
            system('mv traveltimes.tsv S_traveltimes.tsv');   
        end
        
    case 8
        database=varargin{1};
        model_file=varargin{2};
        P_pars_file=varargin{3};
        ratio_pars_file=varargin{4};
        damping=varargin{5};
        relative_weight=varargin{6};
        
        if exist(varargin{7},'file')~=2
            subset=varargin{7};
            locations=varargin{8};
        else
            locations=varargin{7};
            subset=varargin{8};
        end
        
        db = dbopen(database,'r');
        db_site = dblookup(db,'','site','','');
        db_affil = dblookup(db,'','affiliation','','');
        db_arr = dblookup(db,'','arrival','','');
        db_assoc = dblookup(db,'','assoc','','');
        db_orig = dblookup(db,'','origin','','');
        db_event = dblookup(db,'','event','','');

        db = dbjoin(db_orig,db_event);
        
        db = dbsubset(db,'orid==prefor');
        
        db = dbjoin(db,db_assoc);
        db = dbjoin(db,db_arr);
        db = dbjoin(db,db_affil);
        db = dbjoin(db,db_site);

        db_P = dbsubset(db,'phase=~/P/');
        db_S = dbsubset(db,'phase=~/S/');
        
        db_P = dbsubset(db_P,subset);
        db_S = dbsubset(db_S,subset);
        
        ray_make_traveltime(db_P,locations);
        system('mv traveltimes.tsv P_traveltimes.tsv');
            
        ray_make_traveltime(db_S,locations);
        system('mv traveltimes.tsv S_traveltimes.tsv');
        
    otherwise
        help ray_invert_db
        return;
end

runstring = ['raytrace3d PSinvert ' model_file ' ' P_pars_file ' ' ratio_pars_file ' ./P_traveltimes.tsv ./S_traveltimes.tsv 5 5 ' num2str(damping)...
    ' ' num2str(relative_weight) ' 0 180 100 0 360 181 ./inversion_model.model >& invert_output.log'];
system(runstring);
